//
// Created by edward on 22-11-17.
//

#include "miller_rabin.h"
#include "utils.h"

namespace edward {


ll mult(ll x,ll y,ll p)
{
    long double d=1; d=d*x/p*y;
    return ((x*y-((ll)d)*p)%p+p)%p;
}

ll quick_pow(ll a,ll b,ll p)
{
    a%=p; ll ret=1;
    while(b)
    {
        if(b&1) ret=mult(ret,a,p);
        a=mult(a,a,p); b>>=1;
    }
    return ret;
}

bool Miller_Rabin(ll n)
{
    if(n<2) return false;
    const int times=14;
    const ll prime[times]={2,3,5,7,11,13,17,41,61, 73,  101, 97, 7861, 99839};
    for(int i=0;i<times;i++)
        if(n==prime[i]) return true;
        else if(n%prime[i]==0) return false;
    ll x=n-1; while(~x&1) x>>=1;
    for(int i=0;i<times;i++)
    {
        ll a=prime[i]; ll now=quick_pow(a,x,n); ll last;
        if(x==n-1)
        {
            if(now!=1) return false;
        }
        else
        {
            while(x!=n-1)
            {
                x<<=1; last=now; now=mult(now,now,n);
                if(now==1)
                {
                    if(last!=1 && last!=n-1) return false;
                    break;
                }
            }
        }
    }
    return true;
}
}

namespace gfg {
    //https://www.geeksforgeeks.org/primality-test-set-3-miller-rabin/
    // C++ program Miller-Rabin primality test
    using ll = long long;

// Utility function to do modular exponentiation.
// It returns (x^y) % p
    ll power(ll x, unsigned int y, ll p)
    {
        ll res = 1;	 // Initialize result
        x = x % p; // Update x if it is more than or
        // equal to p
        while (y > 0)
        {
            // If y is odd, multiply x with result
            if (y & 1)
                res = (res*x) % p;

            // y must be even now
            y = y>>1; // y = y/2
            x = (x*x) % p;
        }
        return res;
    }

// This function is called for all k trials. It returns
// false if n is composite and returns true if n is
// probably prime.
// d is an odd number such that d*2 = n-1
// for some r >= 1
    bool miillerTest(ll d, ll n)
    {
        // Pick a random number in [2..n-2]
        // Corner cases make sure that n > 4
        ll a = 2 + edward::Random::rand(n - 4);

        // Compute a^d % n
        ll x = power(a, d, n);

        if (x == 1 || x == n-1)
            return true;

        // Keep squaring x while one of the following doesn't
        // happen
        // (i) d does not reach n-1
        // (ii) (x^2) % n is not 1
        // (iii) (x^2) % n is not n-1
        while (d != n-1)
        {
            x = (x * x) % n;
            d *= 2;

            if (x == 1)	 return false;
            if (x == n-1) return true;
        }

        // Return composite
        return false;
    }

// It returns false if n is composite and returns true if n
// is probably prime. k is an input parameter that determines
// accuracy level. Higher value of k indicates more accuracy.
    bool isPrime(int n, int k)
    {
        // Corner cases
        if (n <= 1 || n == 4) return false;
        if (n <= 3) return true;

        // Find r such that n = 2^d * r + 1 for some r >= 1
        int d = n - 1;
        while (d % 2 == 0)
            d /= 2;

        // Iterate given number of 'k' times
        for (int i = 0; i < k; i++)
            if (!miillerTest(d, n))
                return false;

        return true;
    }

}
